/* benchmark guru interfaces */
#include <Configuration.h>
#include <mpi/communicator.h>
#include <mpi/collectives.h>
#include <Utilities/Timer.h>
#include <OhmmsPETE/OhmmsMatrix.h>
#include <Utilities/RandomGenerator.h>
#include <io/dataview.h>
using namespace std;
using namespace APPNAMESPACE;

int main(int argc, char** argv)
{
  mpi::environment env(argc,argv);
  mpi::communicator mycomm;
  OhmmsInfo ohmms("dnsio",mycomm.rank(),0,1);

  //dimensions
  const int D=3;
  ///number of iterations
  int niters=10;
  int gdims[D], parts[D];

  for(int i=0;i<D;++i) gdims[i]=512;
  for(int i=0;i<D;++i) parts[i]=1;
  int npx=8;
  int ic=0;
  while(ic<argc)
  {
    std::string a(argv[ic]);
    if(a == "opt_s") 
      for(int i=0; i<D;++i) gdims[i]=atoi(argv[++ic]);
    else if(a == "opt_p")
      npx=atoi(argv[++ic]);
    else if(a == "opt_i")
      niters=atoi(argv[++ic]);
    //else if(a == "opt_debug")
    //  debug=true;
    ++ic;
  }
  app_log() << " DIMS = " << gdims[0] << " " << gdims[1] << " " << gdims[2] << endl;

  const bool iscomplex=false;
  const int nvars=3;
  typedef float value_type;

  int col_id=mycomm.rank()%npx;
  int row_id=mycomm.rank()/npx;

  {
    mpi::communicator y_slab=mycomm.split(row_id,col_id);

    int npy=mycomm.size()/npx;
    int myslab[D];
    myslab[0]=gdims[0];
    myslab[1]=gdims[1]/npy;
    myslab[2]=gdims[2];

    //create my view
    DataView<3> myview(y_slab);
    myview.init(myslab,npx);
    //for(int i=0; i<mycomm.size(); ++i)
    //{
    //  if(i == mycomm.rank()) 
    //  {
    //    cout << "global rank " << mycomm.rank() << " ";
    //    myview.print(cout);
    //  }
    //  mycomm.barrier();
    //}

    Matrix<value_type> A(nvars,myview.size());
    A=mycomm.rank();
    char fname[16];
    sprintf(fname,"yslab.c%d",row_id);
    Timer myclock;
    double dt=0.0;
    for(int iter=0;iter<niters;++iter)
    {
      myclock.restart();
      write(fname,myview,A.data(),nvars,false,true);
      dt+=myclock.elapsed();
      A=mycomm.rank();
    }
    dt/=static_cast<double>(niters);

    //if(y_slab.rank()==0)
    //{
    //  cout << "global rank " << mycomm.rank() << " row_id " << row_id << " " << dt << endl;
    //}
    double dt_tot=0.0;
    mpi::reduce(mycomm,dt,dt_tot);
    dt=dt_tot/static_cast<double>(mycomm.size());

    double ntot=static_cast<double>(mycomm.size()*A.size()*sizeof(value_type))/1024./1024./1024.;
    app_log() << "#procs #files Size(GB) rate(GB/t) CPU(sec) " <<endl;
    app_log() << "SLABY " << mycomm.size() << " " << npy << " " << ntot << " " << ntot/dt<< " " << dt<< endl;
  }

  {
    mpi::communicator x_slab=mycomm.split(col_id,row_id);

    int npy=mycomm.size()/npx;
    int myslab[D],myparts[D];
    myslab[0]=gdims[0]/npx;
    myslab[1]=gdims[1];
    myslab[2]=gdims[2];
    myparts[0]=1;
    myparts[1]=npy;
    myparts[2]=1;

    //create my view
    DataView<3> myview(x_slab);
    myview.init(myslab,myparts,1);
    //for(int i=0; i<mycomm.size(); ++i)
    //{
    //  if(i == mycomm.rank()) 
    //  {
    //    cout << "global rank " << mycomm.rank() << " ";
    //    myview.print(cout);
    //  }
    //  mycomm.barrier();
    //}

    Matrix<value_type> A(nvars,myview.size());
    A=mycomm.rank();
    char fname[16];
    sprintf(fname,"xslab.c%d",col_id);
    Timer myclock;
    double dt=0.0;
    for(int iter=0;iter<niters;++iter)
    {
      myclock.restart();
      write(fname,myview,A.data(),nvars,false,true);
      dt+=myclock.elapsed();
      A=mycomm.rank();
    }
    dt/=static_cast<double>(niters);

    //if(y_slab.rank()==0)
    //{
    //  cout << "global rank " << mycomm.rank() << " row_id " << row_id << " " << dt << endl;
    //}
    double dt_tot=0.0;
    mpi::reduce(mycomm,dt,dt_tot);
    dt=dt_tot/static_cast<double>(mycomm.size());

    double ntot=static_cast<double>(mycomm.size()*A.size()*sizeof(value_type))/1024./1024./1024.;
    app_log() << "#procs #files Size(GB) rate(GB/t) CPU(sec) " <<endl;
    app_log() << "SLABX " << mycomm.size() << " " << npx << " " << ntot << " " << ntot/dt<< " " << dt<< endl;
  }
  ////myview.init(sizeN,iscomplex);

  //int ng_data=locN*locM*sizeN;

  //Matrix<value_type> A(nvars,ng_data);
  //app_log() << "BW sizeN= " << sizeN << " px= " << nx << " py= " << ny << endl;
  //app_log() << "BW LocN= " << locN << " LocM= " << locM << " localsize= " << A.size() <<endl;
  //A=mycomm.rank();
  ////for(int i=0; i<A.size(); ++i) A(i)=Random();

  //double nbytes=static_cast<double>(A.size()*mycomm.size()*sizeof(value_type))/1024./1024./1024.;
  //hsize_t buffsize=sizeN;

  //Timer myclock;
  //write("cube.0",myview,A.data(),nvars,false,true);
  ////write_pencils_h5(A.cols(),nvars,buffsize,A.data(),"fun",mycomm);
  //double dt=myclock.elapsed();
  //app_log() << "BW col+nobuff GBtype=" << nbytes <<  " time= " << dt << " rate=" << nbytes/dt << " GBytes/sec" << endl;

  //myclock.restart();
  //write("cube.1",myview,A.data(),nvars,true,true);
  //dt=myclock.elapsed();
  //app_log() << "BW col+buff   GBtype=" << nbytes <<  " time= " << dt << " rate=" << nbytes/dt << " GBytes/sec" << endl;

  //char fname[32];
  //myclock.restart();
  //for(int i=0; i<nvars; ++i)
  //{
  //  sprintf(fname,"cube.v%i",i);
  //  write(fname,myview,A[i],1,true,true);
  //}
  //dt=myclock.elapsed();
  //app_log() << "BW col+buff(6) GBtype=" << nbytes <<  " time= " << dt << " rate=" << nbytes/dt << " GBytes/sec" << endl;

  ////myclock.restart();
  ////write_pencils_h5(A.cols(),nvars,A.data(),"fun",mycomm);
  ////dt=myclock.elapsed();

  ////myclock.restart();
  ////write_pencils_h5(A.size(),1,A.data(),"bigfun",mycomm);
  ////dt1=myclock.elapsed();
  ////app_log() << "BW H5 GBtype=" << nbytes <<  " time= " << dt << " rate=" << nbytes/dt << " GBytes/sec" << endl;
  ////app_log() << "BW H5 GBtype=" << nbytes <<  " time= " << dt1 << " rate=" << nbytes/dt1 << " GBytes/sec" << endl;

  return 0;
}

